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ABSTRACT 

Context. The depletion of heavy elements in cold cores of interstellar molecular clouds can lead to a situation where deuterated forms 
of Hj are the most useful spectroscopic probes of the physical conditions. 

Aims. The aim is to predict the observability of the rotational lines of H 2 D + and D 2 H + from prestellar cores. 

Methods. Recently derived rate coefficients for the + H2 isotopic system were applied to the "complete depletion" reaction scheme 
to calculate abundance profiles in hydrostatic core models. The ground-state lines of H 2 D + (o) (372 GHz) and D 2 H + (p) (692 GHz) 
arising from these cores were simulated. The excitation of the rotational levels of these molecules was approximated by using the 
state-to-state coefficients for collisions with H 2 . We also predicted line profiles from cores with a power-law density distribution 
advocated in some previous studies. 

Results. The new rate coefficients introduce some changes to the complete depletion model, but do not alter the general tendencies. 
One of the modifications with respect to the previous results is the increase of the abundance at the cost of other isotopologues. 
Furthermore, the present model predicts a lower H 2 D + (o/p) ratio, and a slightly higher D 2 H + (p/o) ratio in very cold, dense cores, 
as compared with previous modelling results. These nuclear spin ratios affect the detectability of the submm lines of H 2 D + (o) and 
D 2 H + (p). The previously detected H 2 D + and D 2 H + lines towards the core I16293E, and the H 2 D + line observed towards Oph D can 
be reproduced using the present excitation model and the physical models suggested in the original papers. 

Key words. ISM: clouds - Radio lines: ISM - Astrochemistry - Radiative transfer - ISM: abundances - ISM: molecules 
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The rate coefficients for reactive and inelastic collisions between 
the H3 ion and H2 with all possible deuterated variants and 
nuclear spin symmetries have been recentl y derived an d a few 
of them have been experimentally tested bv lHugo et alJ (2009). 
Besides making state-specific astrochemical models viable, the 
cross-sections derived in this work make it possible to model the 
intensities of the dipole-allowed rotational transitions of H_2D + 
and D 2 H + . 

The H2D + and D2H + ions are potentially useful probes of 
pre-protostellar cores (see e. g. Casell Tet al] l2003t IVastel et"aTI 
■ 12004 Ivan der Tak et alJ2 005: Pa gani et alj2009l) . Together with 
the isotopologues and D^, they are likely to belong to the 
most abundant ions in the dense and cold nuclei of prestel- 
lar cores, where the usual tracer molec ules can have prac- 
tically disappeared f rom the gas phase dRoberts et alJ 2003; 
Walmsley etal"1l2004l) - this state of matters is called "complete 
depletion" in the latter paper. While it has been observed that 
heavier substances such as CN ca n exist in the gas phas e at den- 
sities of the order of 10 6 cm -3 dHilv-Blant et al.ll2008h . the as- 
sumption of complete depletion remains a valid first approxima- 
tion at high densities. 

The astrophysical significance of these ions is related to the 
facts that H^ originates almost directly in the cosmic ray ion- 
ization of H2, a nd that it initiates the ion- molecule reactions 
in dense clouds dHerbst & Klemperenll973l) . Furthermore, the 
deuteration of H^ depends on the ortho-para ratio of H2, which 
is likely to be non-thermal and evolve with time in interstel- 



lar clouds dPagani et al.lll992t iGerlich et al] 12001 iFlower et al] 
120061; iPagani et alj|2009l) ~ 

In this paper we apply the newly derived chemical rate co- 
efficients to the "comp letely depleted" case first discussed by 
IWalmslev et al] d2004 . and use the state-to-state coefficients 
to predict the ground-state rotational lines of ortho-H2D + and 
para-D2H + from prestellar cores. The two core models used in 
these simulations correspond to the observed properties of the 
co res Oph D a nd I1 6293E in Ophiuchus. The rate coefficients 
of Hugo et al. (2009)) for so me important deuter ation reactions 
differ from those adopted in lRoberts et alJ (|2003), and in the se- 
ries of papers by Flower, Pineau des Forets & Walmsley, and we 
discuss changes i mplied to the complete depletion model. The 
i Hugo et al] 120 09) rate coefficients have been previously used by 
iPagani et al] ([2009) in the modelling of the prestellar core L183, 
where also CO, N2 and their derivatives were included in the 
chemical reaction scheme. The organization of the present paper 
is as follows: in Sect. [2 we describe the physical and chemical 
models and, in Sect. [3] we present the modelling results. The re- 
sults are discussed in Sect. [4] and finally, the main points of this 
discussion are summed up in Sect. [5] 

2. Model 

In what follows, we describe the model assumptions. These in- 
clude the physical core model, the description of the dust grain 
component, the chemical reaction scheme, and the adopted rate 
coefficients. We also give a brief account of the methods and pro- 
grams used in solving the chemical abundances and the emitted 
molecular line radiation. 



2 



O. Sipila et al.: Modelling line emission of deuterated Hg from prestellar cores 



2.1. Core model 

For the core model we used a modified Bonnor-Ebert sphere 
(BES), which is a non-isothermal clo ud in hydrostat ic equilib- 
rium (previously discussed by, e.g.. Evans et alJuOOlh . heated 
externally by the interstellar radiation field, ISRF. The possi- 
bility of a hydrostatic or near-equilibrium configuration is sug- 
gested by observations towards some prestellar cores which 
seem to represent an advanced stage of chemical evolution (char- 
acterized by a high degree of molecular depletion) , and have 
near-thermal linewidths (e.g., Bergin & Tafalla 20071). 

We first c alculated a den sity profile resulting from an isother- 
mal model (IBonnorl|l956l) . The density pro file was then fed 
into a Monte Carlo radiative transfer program dJuvela & Padoanl 
2003; Juvela 2005) for the temperature calculation. For this pur- 
pose the core model was divi ded into concentric sh ells. We used 
the grain size distribution of lMathis et all d 1977b (MRN) with 
two types of dust grains: carbonaceous and silicate (the grain 
model is further discussed in Sect. \2.2A\ . The choice of using 
the MRN distribution implies that the abundance of very small 
grains (with radii less than ~ 0.01 /zm) is assumed to be negli- 
gible, and that grain coagulation is not considered. The temper- 
ature profile was calculated separately for the two dust types. In 
each position of the core, the final temperature profile was taken 
to be the arithmetic average of the temperatures of both grain 
types. 

The gas and dust temperatures were assumed to be equal. 
While this is probably a good a pproximation in the dense (> 
10 6 cnr 3 ) central part of the core (Bur ke & Hollenbachl 1983b . it 
may not hold as w ell in the outer layers of the core (~ 10 5 cm 3 , 
Bergin et al. 2006). In the present analysis, which concentrates 
on the very dense nucleus of the core, we have ignored this 
phenomenon. The model core is assumed to be embedded in a 
molecular cloud with extinction corresponding to Ay = 10 m . 
This assumptio n agrees with observations from the vicinity of 
Oph D dChapman et al.ll2009l) . 

To determine the exact core density profile using the com- 
puted t emperature d istribution, we slightly modified the method 
used in lBonnori d!956l) by making the following substitutions: 
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where (3 = k^J ' c / 'AnGm, T c is the temperature in the cent er of the 
core a nd m is the mass of the hydrogen molecule. As in lBonnori 
(fl956l) . A was chosen to be the central density. We also imposed 
the boundary conditions ij/ — 0, di///d^ = and dT/d^ = in the 
centre of the core. Substituting the above expressions into the 
equation of hydrostatic equilibrium (Bonnor's equation 2.2) and 
keeping in mind that T = T(r), one obtains 
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Equation (0 was integrated numerically and the resulting den- 
sity profile was fed back into the radiative transfer program to 
calculate a new temperature profile. This procedure was repeated 
a few times until the density profile converged. 

2.2. Chemistry model 

We adopted t he complet e depletion model discussed in 
IWalmslev etail d2004 and iFlower et all d2004l) . The reaction 



scheme is fairly simple: the gas-phase reactants include only H, 
H + , H2, Hj, and their deuterated forms, He, He + , and free 
electrons. Because of the low temperature associated with the 
environment, the zero-point energies become relevant and thus 
the deuterated forms and different nuclear spin modifications 
have to be considered explicitly. The reaction set also includes 
cosmic ray ionization, H2 formation on grains, and some other 
grain processes. We next discuss the gas-phase and grain-surface 
reactions separately. 



2.2.1. Gas phase chemical reactions 

In all our models, we assumed that the core is electrically neu- 
tral, i.e. that electron abundance equals the difference between 
the total abundance of positive ions and the abundance of nega- 
tively charged grains. 

The chemical reactions and the as sociated rate coe f ficient s 
were compi le d usin g four sourc e s: IWalmslev et al.l J2 004). 
Flow er et al] d2004 . iHueo et al.l d2009t). and iPaeani et al.l 
(I2009I) . The reaction set of Walmsle v et al.l d2004l) was com- 



plemented by IFlower et alJ d2004t) to include the nu clear spin 
modifications of D 2 , D 2 H + , and D^. lHugo etaTI (I2009I) took the 
analysis of the + H2 reacting system further by using a mi- 
crocanonical approach to derive st ate-to-state ra te coefficients. 
This approach differs from that of IFlower et al.l d2004t) in that 
all internal states of the reacting system are explicitly taken into 
account, resulting in detailed state-to-state rate coefficients for 
different configurations of the system. We adopted their ground 
state-to-species rate coefficients for all reactions relevant to the 
Ht + H2 system (including the de uterated forms) and su bsti- 



tuted them int o the c ombined data of lWalmslev et al .I d2004l) and 
Flow er et al. (2004). Furthermore, for dissociative recombina- 
tion (DR) reactions of H3 and its deuterated forms we used the 
newly calculated r ate coefficients presented in Appendix B of 
Pagan i et alJd2009l) . 



2.2.2. Grain reactions 

The H 2 and HP molecules are mainly formed on grai n surfaces 
dGould & Salt)eterill963l:lHollenbach & Salneterl 1971b . whereas 
for the formation of D2, gas-phase reactions are more impor- 
tant. Nevertheless, the formation of all three isotopologues on 
dust grains, and the corresponding destruction terms for atomic 
H and D were included in the chemical system. We assumed that 
ortho and para forms are produced according to their statistical 
population ratios 3:1 and 2:1 for H2 and D2, respectively. 

Grain processes, in particular the formation of ortho-H2 and 
the attachment and destruction of positive ions on grain surfaces, 
are important for the overall degree of deuteration of the cloud 
(see Sect.|4|. The details of these processes (grain surface prop- 
erties, desorption mechanisms, etc.) were not considered in the 
present model. It was assumed that chemical reactions on grains 
(along with the associated adsorption and desorption) take place 
instantaneously. 

The rate coefficients for grain reactions (including, e.g., 
electron attachme nt and recombi n ation of positive ions on 
grains), listed in IWalmslev et al] {2004), were taken from 
IFlower & Pineau de s Foretsl d2003b . in which an MRN distribu- 
tion with fl m ; n = 0.01 //m, a max = 0.3 fim was assumed. These 
limits correspond to an effective grain radius of a e ff = 0.0202 
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Fig. 1. Optical thickness per column densit y according to the 
opacity data of Ossenkop f & Hennin 3 (11994. denoted in the fig- 
ure by 'OH 1994', (solid curve) and lLi &Drainel d200ll) . de- 
noted by 'LD 2001' (dashed curve). The modified extinction 
curve (see text) is also plotted, denoted by 'LD mod' (dotted 
curve). 



/irrQ. Since the rate coefficien t for a grain reaction is proportiona l 
to cr g = na^ s , we scaled the lFlower & Pineau des Foretsl (T2003) 
coefficients by (a e s / 0. 02 fim) 2 , where a e s is the effective grain 
radius of the assumed size distribution. The Coulomb factor, 
J(t, v), which takes into account the electric interaction between 
dust grains and gas phase neutrals and/or ions, was included as 
an additio nal factor. Its mathema tical form has be en discussed 
in deta il in iDraine & Sutinl {l987), and recently in Pagan ieTaT] 
(200§. 



2.2.3. Molecular line radiation 

The abundance distributions of chemical species were deter- 
mined by a chemistry program in physical conditions corre- 
sponding to those in different parts of the model core - in prac- 
tice in each of the concentric shells used in the temperature 
calculations. The radial H2D + (o) and D2H + (p) abundance dis- 
tributions, together with the density and temperature profiles 
were used as input for a Monte Carlo radiative transfer program 
dJuvelall 19971) to predict observable line emission. 

The excitation of the rotational transitions of H2D + (o) 
and D2H + (p) in collisions with para and ortho H2 were 
calculated using the state-to-state rate coefficients from 
iFfugo et alJ d2009h . The data concerning the line frequencies an d 
Einstein A coefficients were obtained from Miller et al. (1989), 
iRamanlal & Tennvsonl (120041) . and lAmano & Hiraol (120051) . 



2.2.4. Dust grain model 

The MPvN model was adopted for the grain size distribution. In 
the density and temperature profile calculations, we used two 
types of grains; carbonaceous and silicate. The optical proper- 



1 The rate coefficie nts of grain reactions listed in Table A.l. of 
Walm slev et alJ i2004l) correspond to a eff = 0.02 //m, contrary to what is 
said in the caption. 
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Fig. 2. Steady state abundances of electro ns, protons, Hi and 
its deuterated form s, calculated usin g the iFlower eta 1. (2004) 
(upper panel) and Hug o et al.l d2009l) (lower panel) rate coeffi- 
cients in an isothermal (T = 10K) model with single size grains 
(a — 0.1 /mi). A cosmic ray ionization rate of f = 3 x 10~ 17 s _1 
is assumed. 



ties of these components were adopted from lLi & Drainel(l200ll) . 
The densities of the carbonaceous and silicate grains were taken 
to be 2.5 gcirT 3 and 3.5 gcrrT 3 , respectively. In the chemistry 
model, no separation between the grain materials was done, and 
a grain material density of 3.0 g cm" 3 was assumed. For the com- 
parison with the lFlowere t al. (2004) results (Sect. 3.1) we used, 
however, the same assumptions of the grain mate rial density (2.0 
gem -3 ) and dust-to-gas mass ratio (0.013) as in Wal mslev et al.l 
(2004). 

The model of Li & Draind d2001l) describes dust in diffuse 
medium. In dense clouds, coagulation and growth of icy grain 
mantles are expected to modify the grain properties and, accord- 
ing to the models of Ossenkopf & Henningf d 19941) . the extinc- 
tion curve is expected to become flatter in the far-infrared and at 
longer wavelengths. Since in our cores the central densitie s are 
high, n(H) ~ 10 6 cm -3 , we modified the lLi & Draineld2001l) dust 
opacities accordingly. This leads to temperatures of the order of 
6 K in the core center. The curves are plotted in Fig.Q] 



3. Results 

In this Section we discuss the differences arising from the use 
of the modified chemical reaction rate coefficients used in this 
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Fig. 3. Steady-state abundance ratios of the nuclear spin modi- 
fications of H2, D2, Hi and its deu t erated forms as functions of 
n(H2) calculated using Hugo et al.l d2009t) rate coefficients. The 
model parameters are the same as used in Fig. [2] 



2.5 




/ = 






2.0 




/p-H 2 D + / - 






/ - 


1.5 




/ - 


1.0 




\ 10 4 o/p-H 2 ^/ 1 








0.5 




o/p-H 3 + 


0.0 








T [K] 

Fig. 4. Steady-state o/p ratios of H2 (multiplied by 10 4 ), H3 , and 
H2D + (top panel), and the steady-state D2H + (p/o) and (m/o) 
ratios (bottom panel) as functions of the gas temperature in the 
range 4 - 20 K. The assumed gas density and grain radius are 
n(H2) = 10 6 cirT 3 and a — 0.1 /mi, respectively. 



paper compared to those used in iFlower eta il (12004 . We then 
present results of chemical modelling carried out using a hydro- 
static core model. We conclude the section with a brief discus- 
sion on simulated line emission spectra. 



3.1. Comparison with \Flower et al\ \2004) results using 
homogeneous models 

We ran a se t of homogeneo us models cor responding to thos e 
presented in lWalmslev et al.l (|2004|) and in lFlower et alj d2004). 
The steady state fractional abundances of the principal ions as 
functions of the gas densit y are shown i n Fig. El The upper di- 
agram corresponds to the IFlower et alj (2004) model and is a 
rep roduction o f thei r Fig. 1. The lower plot is produced using 
the iHugo et all (l2009h rate coefficients for t he H^ + H? i sotopi c 
system and the DR rate coefficients from Pagani et al. (2009). 
The figures were produced assuming an isothermal (T = 10 K) 



model with single sized grains (a = 0.1 /mi) and a grain material 
density of 2.0 g cm 4 . 

The abundance ratios of the most important nuclear spin 
variants as functions of the gas density are shown in Fig. [3] 
and the temperature dependence at low temperatures is illus- 
trated in Fig, ffl The fo rmer should be comp ared with Fig. 
2 of Flower et al.l (120041) and with Fig. 3 of IWalmslev et all 
(12004 , an d the latter with Fig s. 5 and 6 in lFlower et alJ d2004 . 
Following IHugo et al.1 (|2009) we have assigned meta-D^ with 
the modification having the lowest ground state energy, corre- 
sponding to the Ai representation of the symmetry group S 3, and 
ortho-Dj with the E representation (see their Sect. II.A.l). The 
ortho and meta appelations for are th e refore interchanged 
with r espect to those used in IFlower et al. J 20041). iPaeani et ail 
(2009), and most other papers. As in IFlower et all d2004 . the 



D3 (o/m) and D2 (p/o) ratios fall rapidly with increasing n(H2). 
However, the D2H + (p/o) ratio is nearly constant unlike in the 
previous study. The HJ (o/p) ratio is about three times lower 
in the present model. There are also a couple of differences in 
the temperature dependence (Fig. 4). The H2D + (o/p) ratio lev- 
els off be tween 2-3 when ap proaching very low temperatures. 
Using the IFlower et al. (2004) rate coefficients the correspond- 
ing value is ~ 10 in the same conditions. Secondly, the D2H + 
(p/o) reaches sl ightly higher (0.3 - 0.4) at very low tempera- 
tures than in the IFlower et al.1 (2004) model. We note that as our 
model follows the chemistry until steady st ate, the H2 (o/p) rati o 
becomes much lower than in the models of Pagani et al. (2009). 
As a consequence, we obtain a clearly higher degree of deuter- 
ation of H^ and somewhat different nuclear spin ratios for its 
isotopologues as compared with this previous study. We return 
to this matter briefly in Sect. 4. 1 . 

All the plots presented so far assume single size d grains with 
a radi us of 0.1 /mi. As discussed in Sect. 3.3 of IFlower et al.l 
(2004), the grain size has a considerable effect on the abun- 
dances and nuclear spin ratios. The D2 (o/p), D2H + (o/p), and 

(m/o) ratios peak near the grain radius a — 0.1 /mi (see Fig. 9 
in lFlower et al. 2004), whereas the H 2 (o/p), H+ (o/p), and H 2 D + 



(o/p) ratios increase monotonically towards smaller grain radii. 



3.2. Hydrostatic models 

We studied the influence of the grain size distribution and the 
cosmic ionization rate on the chemistry of a depleted, hydro- 
static core. The model BE sphere was assumed to have a central 
density of n(H2) = 2 x 10 6 cirT 3 and an outer radius of 2400 
AU, where the density drops to about 10 5 cm 4 . These parame- 
ters correspo nd roughly to the d ense nucleus of the Oph D core 
according to lHariu et alj d2008l) . The effective grain radius was 
varied by changing the lower limit of the grain size distribution 
and keeping the upper limit constant at 0.3 /mi. Modelling was 
carried out with effective grain radii of 0.05 /mi, 0.1 /mi and 0.2 
/urn. The slope of the extinction curve (Fig. [TJ was modified to 
produce a central core temperature of about 6 K (see Sect. 12.2.41 . 
The chemical reaction network was integrated until t — 2 x 10 6 
years, which was well within steady state (all the abundances 
settled into steady state at times < 10 6 years). 

Figure |5] shows the core density profile corresponding to ef- 
fective grain radius a e ff = 0.1 /mi. The temperature profiles for 
three different grain sizes are superposed. Changing a e ff has a 
marked effect on the temperature profile, but alters the density 
distribution only slightly (the central density is kept fixed in our 
models). The core gets colder as the effective grain radius gets 
smaller; this is because the total grain surface area (effectively 
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Fig. 5. The radial core density profile corresponding to effective 
grain radius a e s — 0. 1 fim, with a central temperature of ~ 6.4 K. 
Radial temperature profiles corresponding to different values of 
effective grain radius (indicated in the Figure) are superposed. 



the total opacity) grows larger as the effective grain radius de- 
creases. 

Figure [6] shows the steady state fractional abundance profiles 
of the principal ions for the selected values of effective grain ra- 
dius and for a cosmic ray ionization rate of f = 3 x 10~ 17 s" 1 . 
The grain radius decreases from top to bottom. As shown in Fig. 
[5] this decrease is accompanied with a slight drop in the average 
temperature. The fractional ion abundances generally decrease 
towards the core centre, i.e., towards higher densities, in accor- 
dance with the results from homogeneous models shown in Fig. 
[2] The density dependence of is less marked, which leads 
to an increased degree of deuterium fractionation in the cen- 
tre when the effective grain radius decreases. The fractional H + 
abundance decreases drastically with decreasing grain size. This 
is caused by accentua ted recombinations on negatively charged 
grains as explained in lWalmslev et al.l (|2004|) . 

The radial profiles of the H 2 (o/p), H+ (o/p), H 2 D + (o/p), 
D 2 H + (p/o), and D3 (o/m) ratios are presented Figure [7] While 
most ratios have rather flat distributions, the D3 (o/m) drops 
towards the centre, i.e. towards higher densities. The density 
dependence is, however, smoother for small grain sizes (bot- 
tom panel). The ortho/para ratios of H 2 , H3 , and H 2 D + increase 
from top to bottom, along with the decreasing grain size. This 
is caused by an intensified replenishment of H 2 (o) when the to- 
tal grain surface area becomes larger. The slight increase of H 2 
(o/p) towards the outer parts reflects the decreasing density (see 
Fig EJ. The H3 (o/p) ratio follows the H 2 (o/p) ratio closely. The 
H 2 D + (o/p) ratio correlates with the H 2 (o/p) ratio, but rises to- 
wards low temperatures (see Fig. 0). These tendencies result in 
a slight increase of the H 2 D + (o/p) ratio towards the core centre. 
The D 2 H + (p/o) ratio changes very little as a function of radius 
and from model to model. Like the D3 (o/m) rat io, the D 2 H + 
(o/p) ratio has a flat maximum around a ~ 0.1 /mi dFlower et alj 
2004), and depends weakly on the density (Fig. [3]). 

The effect of changes in the cosmic ionization rate, is il- 
lustrated in Fig. [8] This shows the radial distributions of the 
principal ions for the cases where £ is made 10 times lower 
(top) and 10 times higher (bottom) than in the models discussed 



above. Lowering £ decreases the degree of ionization. An en- 
hanced cosmic ray ionization increases the ion abundances, but 
decreases the deuteration through an intensified e lectron recom- 
bination ( Wal mslev et alJ2004tlFlower et alJ2004l) . The effect is 
particularly marked for D3 near the outer edge of the core. 

3.3. H 2 D + and D 2 H + spectra 

We calculated the spectral line profiles of the ground state tran- 
sitions of H 2 D + (o) (lio - In, 372 GHz) and D 2 H + (p) (1 10 - l i, 
692 GHz) from models corresponding to the observed properties 
of the prestellar cores I16293E and Oph D. 

The very narrow 372 GHz line of H 2 D + (o) detected with 
APEX towards the density maximum of Oph D has been sug- 
gested to originate in a hydrostatic core (Harju et al. 2008). In ac- 
cordance with this suggestion the physical model used for Oph D 
is a hydrostatic core heated externally by the interstellar radi- 
ation field (ISRF). Judging fro m previous maps, in particular 
the ISOCAM 7 fim image from iBacmann et al.l (120001) . we as- 
sumed an angular radius of 20". We adopted the recent distance 
estim ate of 120 pc to the Ophiuchus complex (lLombardi et al.l 
2008), which implies that the adopted radius corresponds to 
2400 AU. The obscuration provided by the surrounding molec- 
ular cloud was assumed to correspond to Ay = 10 m . We varied 
the central H 2 density (no), the effective grain radius (a e ffX an d 
the cosmic ray ionization rate (0, and calculated the appropriate 
abundance profiles using the chemistry model described above. 
Finally, the populations of the rotational levels of H 2 D + (o), and 
the 1 10 - 1 1 1 line profiles were calculated using the state-to-state 
rate coefficient s for H 2 D + (o) + H ? (para and ortho separately) 
collisions from Hu go et all {2009) and our Monte Carlo radia- 
tive transfer program. Calculations with different values for the 
model parameters were carried out until a reasonable agreement 
with the observed H 2 D + (o) profile was met. This model was then 
used to predict the D 2 H + (p) line at 692 GHz as observed with 
APEX. 

The results of our simulations are presented in Fig.[l0] which 
also shows the observed H 2 D + (o) spectrum. The best fit was 
obtained using the following parameters: no — 2 x 10 6 cm~ 3 , 
a e ff = 0.1 fim, and ( = 6x 10~ 17 s" 1 . The model core has a cen- 
tral temperature T = 6.36 K. The fractional H 2 D + (o) abundance 
is ~ 3.6 x 10~ 10 in the core center, rising to ~ 7.2 x 10~ 10 at the 
edge. A Gaussian fit to the model spectrum yields an FWHM of 
0.35 km s -1 ; the peak optical depth is T pea k ~ 1-5. A molecular 
line profile for the D 2 H + (p) 1 10 — loi transition using the same 
core model parameters is shown in the lower panel of the figure. 
A Gaussian fit to the spectrum yields an FWHM of 0.27 km s" 1 , 
the peak optical depth is r pea k ~ 0.6. The D 2 H + (p) abundance is 
~ 1.4 x 10" 10 at the core center and ~ 1.1 x 10" 10 at the edge - 
essentially a flat profile with a minor maximum at around 1200 
AU. 

Both H 2 D + and D 2 H + have been d etected with the CSO to- 
wards I16293E bv IVastel et"ail J2004). To our knowledge this 
work has produ ced the only rathe r firm detection of the 692 
GHz line so far. IVastel et all (120041) estimated the H 2 D + (o) and 
D 2 H + (p) column densities from the observed line profiles by 
assuming line-of-sight homogeneity and an excitation temper- 
ature of T ex = 10 K. We assumed that these column densi- 
ties are approximately valid, and adjusted the fractional abun- 
danc es accordingly. The physical model of the core was adopted 
from Star k et al.l (120041) . allowing for the fact t hat the cloud is 
proba bly nearer than previously thought (120 pc, Lombard i et al.l 
120081) . The model consisted of a compact, homogeneous nucleus 



6 



O. Sipila et al,: Modelling line emission of deuterated Hg from prestellar cores 



with a density n® — 1.6 x 10 6 cirT 3 up to a radius of 750 AU, 
surrounded by an envelope with the density power law n oc r 
extending to a distance of 6000 AU (50"). A dust temperature 
Tdust = 16 K derived by IStark et al.l d2004l) was adopted as the 
temperature of the nucleus. The dust temperature was gradually 
increased toward the core edge, reaching a value T^ust = 20 
K at 60 00 AU. It should b e noted, however, that the observa- 
tions of IVastel et al.l (120041) suggest a lower value for the re- 
gion where the H2D + and D2H + emission originate. Assuming 
thermal broadening only, the line widths listed in their Table 
1 imply kinetic temperatures T^ n = 11.2 + 2.5 (H2D + ) and 
Tkin = 9.1 ± 4.4 (D2H + ). These estimates are upper limits as 
they neglect the possible turbulent broadening and instrumental 
effects. In order t o reproduce the column densities derived by 
IVastel etall d20 04). the fractional H2D + (o) abundance was set to 
~ 10" 10 in the nucleus, and it was let to fall radially to ~ 10~ n at 
the outer edge. The D2H + (p) abundance was set at 0.6 times the 
H2D + (o) abunda nce, which is with in the column density error 
range reported bv lVastel et al] (120041) . 

The line profiles calculated from this model are shown in 
Fig HU The simulated H2D + (o) line has a peak optical depth r 
~ 0.33. A gaussian fit to the spectrum yields an FWHM of ~ 
0.46 kms -1 and a peak antenna temperature T*^ ~ 1.42 K. The 
simulated D 2 H + (p) line yields r ~ 0.30, an FWHM of ~ 0.41 
kms" 1 and a peak T* k ~ 0.44 K. The antenna temperatures of 
the simulated H 2D + (o) and E>2H + ( p) spectra agree well with the 
observations o f IVastel et al.l (|2004), but the lines are too broad, 
probably because the assumed kinetic temperature is too high. 
We return to this issue in Sect. 14.31 



4. Discussion 

We now discuss in more detail the results presented in the last 
Section, starting with the differences in homogeneous models 
arising from using different rate coefficients. We then discuss 
the abundance distributions in our hydrostatic core models, and 
conclude by looking into the simulated line emission spectra. 



4. 1 . Homogeneous models 

In the mod els calculated using the new reaction rate coefficients 
derived by iHugo et al] d2009j). the deuteration o f Hi proceeds 
faster than in the models of iFlower et alJ d2004l) . As a conse- 
quence, becomes the dominant ion at relatively low den- 
sities. The totaQ rate coefficients of IHugo et al.l (I2009I) for the 
deuteration reactions Hi 



H 2 D + -> D 2 H + D+ are larger 



than those used bv IFlower et al .1 d2004l) bv an average factor of 
~ 4-5. This difference origin ates in the discrepa ncy between the 
labora tory measurements of lGerlich et alJ d2002l) and lHugo et alJ 
(2009, see their Sect. IV. C) concerning the forward rate coeffi- 
cients at very low temperatures. On the other hand, the total rate 



2 By the "total rate coefficients" we mean the resultant forward and 
backward rate coefficients, k t and fc_i, in the presentation omitting the 

h 

nuclear spin variants, e.g., Hj + HD ^ H 2 D + + H 2 . In this presentation, 

the H 2 D + formation rate can be written k\ [HJ ][HD], where k t is a com- 
bination of seven individual rate coefficients weighted according to the 
fractions of Hg (o) or Ht (p) of the total Hj abundance (see Table VIII 
of Hugo et al. 2009). The other deuteration reactions and their reverse 
reactions are composed of eight reactions. Besides the temperature, the 
total rate coefficients and the equilibrium constant, K = k/k-, depend 
on the relative populatio ns of nuclear spin variants, in particular o/p H 2 
Ce.g. lGerlich et al j2002h . 
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Fig. 6. Radial fractional (with respect to n(H2)) adundance pro- 
files of electrons, H + , H3 and its deuterated forms at t — 1 x 10 6 
years, for a cosmic ray ionization rate £ : 



coefficients for the "backward" reactions (k-\, fc_2, k^) are sim- 
ilar in the two reaction schemes. This can be traced to the fact 
that the "backward" rate coefficient for the reaction H2D + (o) + 



H 2 (o) ^ 0) H+(o/p) + HD derived bvlHugo et all d2009h agrees 
with the estimate of iGerlich etall (120021 Sec t."3~3). 

The ratios K = k/k- are thus larger in iHugo et al.l (120091) 
which implies that chemical equilibrium is established at higher 
concentrations of deuterated species. The equilibrium constants 
in the present model ar e, ho wever, smaller than in the models 
of iRoberts etail ( 120041) and ICaselli et all d2008l) . The forward 
rate coefficients of iHugo et al. I (12009b are similar to those used 



O. Sipila et al,: Modelling line emission of deuterated from prestellar cores 



7 



tteff = 0.2 fim 



(o/m) 



(P/O) D 2 H+ 

(o/p) H 2 D+~ ~"~ ,L1 ^-u.^ 



1e4 x (o/p) H 2 



(o/p) H+ 



a c ff = 0.1 /tm 



(o/p) H 2 D+ 



1e4 x (o/p) H 




(o/m) D 



(p/o) D 2 H 4 



+ 



+ 



a e g = 0.05 /im 



1e4 x (o/p) H 2 



(o/p) H 2 D+ 



(o/p) H 




(o/m) D 



(p/o) D 2 H 4 



2.5 
2.0 
1.5 
1.0 
0.5 

2.5 
2.0 
1.5 
1.0 
0.5 

2.5 
2.0 
1.5 
1.0 

0.5 



500 



1000 1500 
R[AU] 



2000 2500 
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H2, and its deuterated forms at t — 1 x 10 6 years, for a cosmic 
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in the two previous studies while the backward rate coefficents 
are several orders of magnitude larger owing to the non-thermal 
H2 (o/p) ratio. 

As discussed in lFlower et al.l d2006l) . the ortho/para ratios of 
H3 and H2D + correlate with H2(o/p) through proton exchange 
reactions. The HJ ion is primarily produced in the para form via 
the reaction of HJ (p) with H2(p). The para-ortho and ortho-para 
conversions through reactions with H2(o) are, however, rapid, 
and compete only against the deuteration reaction with HD. In 
the conditions considered here a good approximation of the o/p 
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Fig. 8. Same as the middle panel in Fig. [6] but for cosmic ray 
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ratio can be obtained from 



H 3 + (o/p) 



k po H 2 (o/p) 



k lo x(HD)+k op n 2 {o/p) 



where k po and k op are the rate coefficients of the nuclear spin 
changing reactions of HJ(p) and H3 (o) with H2(o) (reactions 6 
and 7 in Flower 2006), and k\„ is the rate coefficient of the re- 
action of H^(o ) with HP which produces m ainly H2D"*"(o) (see 
Table VIII in iHugo et al]l2009l) . Using the iHugo et al l ([2009) 
rate coefficients, with k po ~ k op ~ 0.3fe[ O , we obtain a steady- 
state ratio of H^(o/p) ~ 0.2 - 0.4 in cold, dense gas, where 
H 2 (o/p)~x(HD)(seeFig.|4j. 

Also for H2D + , the nuclear spin changing reactions with 
H2 are effective and an approximate value of the o/p ratio can 
be obtained t hroug h balancing these according to Eq. (7) of 
iGerlich et al.l (12002). Below ~ 8 K, however, primary production 
from + HD and destruction through the secondary deutera- 
tion reaction H 2 D + + HD -> D 2 H + + H 2 start to dominate. At 
very low temperatures the H2D + (o/p) ratio settles somewhere 
between 2 and 3 (Fig. |4j». 

The production and destruction of ortho and para D2H + oc- 
cur mainly via the deuteration reactions H2D + + H D — » D2H + + 
H 2 an d D 2 H + + HD -» D+ + H 2 . As discussed in iFlower et al.l 
(2006), the lower energy ortho state is favoured in the reverse 
reactions with H2(o) (reaction s 11 and 12 in their Sect. 3.4; 
Table VIII in lHugo et ai1l2009l) . Moreover, D 2 H + (p) is reduced 
through para-ortho conversion in the reaction with HD. In the 
present model the p/o ratio is ~ 0.2 - 0.4 below 10 K and is not 
sensitive to n(H2). 
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We note that while the H2D + (o/p) rati o is similar to that 
found in the models of IPagani et al] ^2009, Fig. 9), the D 2 H + 
(p/o) ratio is larger by a factor of 2-4 in the present model. By 
performing some test runs we could establish that the main rea- 
son for this apparent discrepancy i s in the diffe r ent H2 (o/p) ratio. 
In the time-dependent model of Paga ni et al.l ((2009), the initial 
H2 (o/p) ratio is larger than unity, and the chemical evolution is 
followed until the model matches with the constraints based on 
observations of the N2D + /N2H + ratio toward the L183 core. The 
resulting H 2 (o/p) ratio (~0. 005-0. 05) is much higher than in our 
models (less than 10 4 ). At a high H 2 (o/p) ratio, D 2 H + (p) is 
mainly destroyed in the reaction D2H + (p) + H2(o) — > H2D + (p) 
+ HD, and the destruction by HD (including para-ortho conver- 
sion) is not significant. 

The principal r eactions dete r minin g the ortho/meta ratio of 
D3 are discussed in lFlower et al. (2004) (Appendix B, with inter- 
changed meta and ortho appellations with respect to the present 
work). The (lowest energy) meta form of D3 is primarily formed 
via deuteration of D2H + (o) by HD, and via ortho- meta conver- 
sion in reaction D3 (o) + HD — » D3 (m) + HD. The destruction of 
DJ (m) is overwhelmingly dominated by dissociative recombina- 
tion with electrons or negatively charged grains. For D3 (o), the 
reaction D2H + + HD — > D3 + H2 in the backward direction and 
the ortho-meta conversion mentioned above are comparable with 
electron recombination as de struction mechanism s. The formula 
presented in the appendix of iFlower et al.l (120041) approximates 
the equilibrium o/m ratio with a high accuracy. The ratio is pro- 
portional to the electron density, and decreases as the gas density 
increases. On the other hand, it changes very little as a function 
temperature in the range considered here (see Figs. 3 and 4). 

Finally, we note that the D2 (p/o) ratio correlates with the 
DJ (o/m) ratio (see Fig. 3). The higher energy D2(p) is prin- 
cipally formed from D3 (o) via dissociative recombination with 
free electrons and on dust grains, whereas D2(o) forms primarily 
from D3 (m) (DR rate coefficients from Pagani et al. 2009 and 
branching ratios from Flower et al. 2004 have been used here). 
The destruction of D2(p) is dominated by deuterium exchange 
reactions with D3 (m) and deuteration of H3 (p), while for D2(o) 
the deuteration of D2H + (o) and H^(p) are the most important 
destruction reactions. The resulting equilibrium D2 (p/o) ratio is 
roughly an order of magnitude lower than the D3 (o/m) ratio. 

4.2. Hydrostatic models 

In the hydrostatic models (Figs. [6] and [7} combined effects of 
changing grain size, density and temperature can be seen. The 
decrease of the average grain size from the top panel to the bot- 
tom signifies an increase in the total grain surface area. This re- 
sults in a diminishing fractional abundance of free electrons and 
a growing abundance of negatively charged grains (owing to the 
sticking of electrons onto neutral grains), and an increase of the 
H2 (o/p) ratio (because of enhanced H2 formation on grain sur- 
faces). The radial fall off of the density (from left to right) is 
accompanied by an elevated temperature at the outer edge of the 
core. 

The DJ ion is the most abundant deuterated form of H3 in the 
centre of the core. In the following we use the D3/H3 abundance 
ratio to quantify the degree of deuteration. The D3/H3 ratio 
changes as a function of grain size. For the largest grain size, a = 
0.2 fjm, [D+]/[H+] is ~ 3.3 in the core center, and ~ 0.08 at the 
edge. As the grain size is decreased to a = 0.1 fim, [D^]/[Hi] 
increases both in the center (to ~ 8.6) and at the edge (to ~ 0.6). 
For still smaller grains, a = 0.05 ^m, [D 3 h ]/[H 3 h ] drops slightly 
(to 8.5) in the center, but increases to ~ 1 .4 at the edge. 
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Fig. 9. Dependence of the D fractionation on the average grain 
size. The curves correspond to different values of the gas density. 

The changes in the deuterium fractionation as functions of 
the density and the effective grain radius are plotted in Fig. [9] 
for isothermal models with 7\i n = 7 K. One can see in the plot 
that for each value of density, there is a deuteration maximum 
which shifts to lower values of a as the density is decreased. The 
sloping down of [Dt]/[Ht] on the right, towards larger grains, 
is caused by an increasing abundance of free electrons. The de- 
crease of [DjJ/rHJ] on the left side, towards smaller grains, is 
caused by the increase in the abundance of H2(o) which hinders 
deuteration by converting H2D + back to H3 , and by the increase 
of negatively charged grains which replace electrons as the prin- 
cipal destructors of Dt . 

Figure [8] shows radial abundance profiles in a chemical 
model otherwise similar to the middle panel of Fig. [6] but for 
cosmic ionization rates £ = 3 x 10~ 18 s _1 and f = 3 x 10 16 
s . It can be immediately seen in Fig. [8] as compared with 
the middle panel in Fig. [6] that the abundances of H + and H3 
follow the change in the cosmic ray ionization rate. The H + 
abundance is proportional to £ via dissociation of H2 by cos- 
mic rays and via diss ociative charge transfe r reaction between 
He + and H 2 (see, e.g. lWalmslev et alJl2004 Appendix A). The 
formation of H3 depends on £ through Hj , which is the princi- 
pal produ ct of the collision between H2 and a co smic ray pro- 
ton (e.g., IWalmslev et alJl2004t IPagani et alJl2009t) . The H 2 D + 
and D2H + abundances are dragged upwards with H3 with an in- 
tensified cosmic ray ionization, whereas the D3 abundance re- 
mains roughly constant in the core center and falls off in the 
outer regions. So, quite reasonably, the [Dj]/[HJ] ratio which 
depends on DR reactions diminishes strongly with an increasing 
£ whereas the effects on [H 2 D + ]/[H+] and [D 2 H + ]/[H+] which 
depend on HD and H2(o/p), are less notable. 

4.3. Molecular line profiles 

The observed H2D + (o) profile from Oph D can be reproduced 
reasonably well using a hydrostatic core model together with 
the "co mplete depletion" c hemistry model with the rate coeffi- 
cients of Hugo etafl (120091) (see Fig. ITOb . The model giving the 
best agreement with the observed H2D + (o) spectrum (see Sect. 
133} predicts a very weak D 2 H + (p) line at 692 GHz (7/* ~ 0.07 
K with APEX). This line has a lower peak opacity (r ~ 0.6) 
than the H2D + (o) line, whereas the radial T ex distributions of 
the two lines are similar. One should note that even an optically 
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the D2H + 1 io - loi (para) transitions as observed with APEX for 
the best-fit core model ( solid curves) . Also plotted is the obser- 
vation toward Oph D bv lHariu et ail d2008l) for the H 2 D + 1 10 - 
In transition (dotted curve). 



In) and the D2H + (lio - 1 oi) transitions for the prestellar core 
16293E as observ ed with the CSO, using the core parameters of 
Stark et al.l d2004l) (solid curve). The modified model (see text) 
is also plotted (dashed curve). 



thick D2H + (p) (lio - loi) line from a cold source is weak in the 
brightness temperature scale because the frequency lies far away 
from the range where the Rayleigh- Jeans approximation is valid. 

The physical model adopted from Stark et al. (2004) (see 
Sect. 13.31 1 for the core I16293E pro duces easily the b right 
H 2 D + (o) andD 2 H + (p) lines observed bv lVastel et al.l (12004 (see 
Fig. QT| solid c urves). The relative ly high temperature (16 - 20 
K) assumed by IStark et akl ([2004) gives, however, rise to too 
broad lines as compared with the observations. Another prob- 
lem with this assumption is that the chemistry model predicts a 
clearly highe r D2H + (p)/H2D + (o) abundance ratio around 16 K 
(see Fig. 7 in lFlower et a l. 2004). We carried out another simu- 
lation assuming that the kinetic temperature is 10 K in the center 
and that it increases gradually to 12 K at the edge. The fractional 
U_2D + (o) and D2H + (p) abundances were assumed to be roughly 
equal, in agreement with the prediction of the chemistry model. 
The modelled spectra are shown in Fig.QT](dashed curves). Now 
both intensi t ies and line widths are within the range reported by 
IVastel et all (|2004|) . The H 2 D + (o) and D 2 H + (p) column densi- 
ties, bot h ~ 2.3 X 10~ 13 cm -2 , are slightly higher than those de- 
rived by IVastelltaD d20Q3>. 



5. Conclusions 

We studied the chemistry in molecular cloud cores in the com- 
plete depletion limit, taking advantage of new state-to-state 



chem ical reaction rate coefficients for the H3 + H2 reacting sys- 
tem dHugo et al.ll2009h . We used the modelling results to pre- 
dict line emission spectra for the H2D + 1 10 - 1 1 1 (372 GHz) and 
the D2H + 1 10 - loi (692 GHz) transitions in a core model corre- 
sponding to the prestellar core Oph D. We also carried out sim- 
ulations of the above transitions for the I16293E core . The sim- 
ulated spe ctra were compared with observations by Harju et al. 
(2008) and lVastel et al.l d2004l) for the respective cores. 

In this paper, we used steady-state chemical abundances in 
the molecular line simulations. It has been noted in the literature 
that the timescale for chemical development can be an order of 
magnitude larger than the dynamical timescale, and thus chemi- 
cal steady-s tate might not be rea ched before collapse begins (see 
for example Fl ower et al.ll2006l) . In general, however, ions reach 
chemical steady-state earlier than neutral species, due to the 
speed of the ion-molecule reactions. In fact, our chemical mod- 
els indicate that ions reach steady-state before t — 1 x 10 5 years. 
On the other hand, the relatively high H2D + (o) and D2H + (p) col- 
umn densities derived toward Oph D and I16293E imply an ap- 
preciable degree of deuteration, which suggests that the cores 
are in an advanced chemical state. Taking into account these ar- 
guments, the assumption of steady-state seems, in the present 
context, appropriate. 

The simulated profile of the H2D + (lio - In) line from 
a self-consistent hydrostatic core model (including chemistry 
modelling) is in reasonably good agreement with observations 
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toward Oph D. The model predicts that the D2H + ( 1 10 — loi)line 
at 692 GHz probably cannot be detected in this source. In fact, 
this transition is difficult to detect in emission in any cold object 
because of its high frequency. The same is true for the H2D + (p) 
and D2H + (o) ground state lines, which lie at even higher fre- 
quencies of ~ 1.4THz. 

The Oph D core seems particularly appropriate for using a 
hydrostatic model. However, the chemistry model used here be- 
gins to lose its validity in the outer parts, where the density drops 
below 10 6 cirT 3 . In these conditions one can no longer ignore 
heavier substances in the gas phase, which modify the abun- 
dances of HI and its deuterated forms. So, in reality, the H2D + 
and D2H + abundances which in our models rise toward the core 
edge are likely to turn down outsi de the dense nucle us. This ef- 
fect can be seen in the models of Pagani et al. (2009). It is obvi- 
ous that chemical modelling should be extended to account for 
the presence of heavier species. This would also justify studies 
of larger cores, presenting larger density gradients. 

Our simulations of molecular line profiles toward the 
I16293E core are in good agr e emen t with the observational re- 
sults presented in IVastel et alj d2004l) . While line e mission sim- 
ulation s utilizing the physical core parameters of Star k et ail 
(2004) produce enough emission , we found that to produce 
linewidths comparable to those in IVastel etall (|2004), the core 
temperature profile and chemical abundances should be modi- 
fied. This modification is also justified by results from chemi- 
cal modelling. The H2D + (o) and D2H + (p) column densities are 
found to become comparable when the temperature drops to ~ 
10 K. 

The deuteration of Hi proceeds mu c h fast er in the present 
reacti on scheme using the iHugo et aTl (120091) results than in 
the iFlower et alj ([2004) mo del where the deute ration rate co- 
efficients were adopted from iGerlich et al] d2002l) . As a conse- 
quence, the steady-state D3 /Ht abundance ratio becomes larger 
than in the models of Flower, Pineau des Forets & Wal msley. 
Howev er, it should be noted that already in Wal msley et alj 
(120041) . becomes the most abundant isotopologue of H3 at 
very high densities, and that th is phenomenon was predicted pre- 
viously by iRoberts et al] d2003l) . 

Th e low-temperature rate coefficients from IGerlich et al] 
(2002) are lower than those bv lHugo etall d2009t) by a factor of 
~ 4. The experimenta l circumstances have been thoroughly dis- 
cussed in Hug o et aTl d2009h . concluding that further experimen- 
tal investigations, preferably with a different setup, are urgently 
needed. The ortho/para ratios of the deuterated forms of H3 are 
also modified in the new model. Along with the changes in abun- 
dances, this has an impact on the observability of H2D + and 
D2H + toward prestellar cores. We find at very low temperatures 
a lower H2D + (o/p) ratio, and a slightly higher D? H + (p/o) ratio 
than p redicted by using the coefficients adopted in lFlower et al] 
(2004). 

Th e state-to-state coefficients calculated by IHugo et aTl 
( 2009) provide an opportunity to refine both the chemistry model 
and the model for the excitation of the rotational transitions of 
interest. Besides collisions with para and ortho H2, it would seem 
reasonable to examine the excitation of rotational levels of H2D + 
and D2H + in collisions with HD. The next step forward is to 
combine the full state-to-state reaction network with radiative 
transfer calculations. The realization of this improvement looks 
particularly feasible for the complete depletion model. 
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